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process,  step  by  step,  reflecting  the  research  nature  of  the  program.  On  the 
other  hand,  the  ITERATE  command  allows  the  iterations  to  be  automatic.  Once 
a  STOP  command  is  given,  the  run  will  be  terminated. 

Printed  output  generated  by  FEARS  consists  of  summary  data  after  each 
iteration  step.  The  user  also  has  the  option  to  reprint  initial  input  and/or 
detailed  information  about  the  solution  process  during  the  iteration.  The 
formats  of  these  printouts  are  described  in  Chapter  4.  In  addition,  the  DUMP 
command  generates  a  dump-file  of  the  data  structure  which  can  be  used  as  input, 
either  for  FEARS  itself,  or  for  various  post-processors.  The  format  of  this  dump- 
file  is  described  in  Section  4.7. 

Chapter  5  gives  the  computer  dependent  control  statements  necessary 
to  run  FEARS  on  the  UNIVAC  1100  computer.  The  program  is  supplied  with  a 
dummy  subroutine  giving  zero  values  for  the  bilinear  matrices  E  and  D  , 

(see  Section  2.7).  If  these  matrices  are  coordinate  dependent,  the  user  must 
supply  the  appropriate  subroutine  replacing  the  supplied  dummy  routine.  This 
preparation  must  be  performed  prior  to  the  use  of  FEARS. 
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0.  Introduction 


The  FEARS  program  normally  takes  the  problem's  descriptive  data  from 
the  run-stream.  This  input  data  is  composed  of  the  Geometry  ^Chapier*  1) — ^ — 
and  of  the  Bilinear /Error  matrices^Ghaptcr  Alternatively,  this  data 

3-  ^ 

can  be  read  from  ile  ,  which  has  been  generated  by  a  previous 

execution  of  FEARS  1  (SclIIou  3»  13^f  Once  the  initial  data  is  read,  FEARS 
enters  into  a  command-mode.  In  the  command-mode,  the  user  must  input  a 

command  (with  some  parameters)  instructing  the  program  about  the  next  step 

^ - -  -  .  . .  '  "  - ^ 

to  be  taken.  The  available  commands  are  given  in  Chapter  3.  This  mode  allows 

the  user  to  have  complete  control  over  the  iteration  process,  step  by  step, 
reflecting  the  research  nature  of  the  program.  On  the  other  hand,  the 
ITERATE  command  allows  the  iterations  to  be  automatic.  Once  a  STOP 
command  is  given,  the  run  will  be  terminated. 

Printed  output  generated  by  FEARS  consists  of  summary  data  after 
each  Iteration  step.  The  user  also  has  the  option  to  reprint  initial 
input  and/or  detailed  information  about  the  solution  process  during  the 
iteration.  The  fozmats  of  these  printouts  are  described  in  Chapter  4. 

In  addition,  the  DUMP  command  generates  a  dump-file  of  the  data  structure 
which  can  be  used  as  input,  either  for  FEARS  Itself,  or  for  various  post¬ 
processors.  ^The  format  of  this  dump-file  is  described  in  Section  4.7. 
c  Chapter  5  gives  the  computer  dependent  control  statements  necessary 
to  run  FEARS  on  the  UNIVAC  1100  computer.  The  program  is  supplied  with  a 
dummy  subroutine  giving  zero  values  for  the  bilinear  matrices  E  and  D, 

"faee — Sectidn~2 . 7)i-  If  these  matrices  are  coordinate  dependent,  the  user 
must  supply  the  appropriate  subroutine  replacing  the  supplied  dummy  routine. 
This  preparation  must  be  performed  prior  to  the  use  of  FEARS. 

1  'T 

\ 
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Chapter  I.  Geometry  Input 
1.1  Introduction  —  0-D,  1-D,  and  2-D  domains 

The  domain  V,  on  which  the  problem  is  to  be  solved,  must  be  initially 
broken  up  into  subdomains.  Each  subdomain  is  a  generalized  quadrilateral, 
having  4  corner  points  and  4  sides,  with  each  side  being  either  a  straight 
line  or  an  arc  of  a  circle.  Since  a  unit  quare  will  be  mapped  onto  each 
of  these  subdomains,  they  should  not  be  "degenerate"  or  "singular".  For 
example,  the  angles  formed  at  the  vertices  should  not  be  too  near  0°  or 
180°,  the  domains  should  not  be  nearly  triangular,  and  no  overlapping  of 
the  sides  is  allowed.  Thus,  it  is  good  practice  to  avoid  subdomains  as 
shown  in  Figure  1.1. 

X 

Figure  1.1.  Subdomains  to  be  avoided. 

Figure  1.2  illustrates  how  a  disk,  and  triangular  and  pentagonal 
domains  may  be  substructured  into  generalized  quadrilaterals. 


Figure  1.2.  Subdomain  structuring  for  a  disk,  triangular 
domain,  and  pentagonal  domain. 


Because  FEARS  uses  blending  functions  to  map  the  unit  square  onto 
each  generalized  quadrilateral  (see  FEARS  Mathematical  Description, 

Chapter  I) ,  domains  having  boundaries  composed  of  circular  arcs  are 
represented  exactly  and  are  not  approximated  by  isoparametric  elements. 
Therefore,  the  approximations  of  the  geometry  and  the  exact  solutions 
are  not  mixed,  when  using  FEARS. 

The  corner  points  of  the  subdomains  are  simply  called  points  or 
0-D  domains  (zero  dimensional  domains).  These  are  denoted  by 

k-1.2....N0. 

The  open  line  segments  ioining  the  points  are  called  lines  or  1-D  domains, 
and  are  denoted 

P*  k«l,2,...Nl. 

Finally,  the  open  subdomains,  each  with  4  lines  and  4  points  forming  its 
boundary,  are  called  2-D  domains,  and  are  denoted 

V £  k-1,2, . . ,N2. 


Figure  1.3.  The  substructuring  and  labelm*  of  -  Disk. 
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In  the  geometry  input  the  following  information  must  be  supplied 


(i)  the  number  of  points ,  lines,  and  2-D  domains, 

(ii)  the  global  coordinates  of  the  points, 

(iii)  the  end  point  indices  for  each  line, 

(iv)  the  radius  of  each  line, 

(v)  the  boundary  conditions  for  each  line  and  point, 
and  (vi)  the  comerpoint  indices  of  each  2-D  domain. 

1.2  Format 

The  format  for  the  geometry  input  is  as  follows: 


1.3.  Description  of  the  Format 


NO  Is  the  number  of  points  (O-D's).  The  data,  i,x^,y^,b^,  \ 
describe  the  points  as  follows: 


1  -  the  Index 


of  P°,  1-1,..., NO, 


*  t*ie  x_y  coordinates  of 

b^  -  the  boundary  condition  at  where 


0  *>  u^  free,  U2  free, 

1  ■>  u^  fixed,  u2  free, 

2  ->  u ^  free,  u^  fixed. 


3  *>  u^  fixed,  u2  fixed  J  ,  and 


^Vi*Wi^  *  the  solution  value  (u^.u^)  at  P°.  These  must  be  given 
even  if  the  boundary  condition  is  free,  in  which  case 
any  value  for  v^Wj  may  be  given. 

N1  is  the  number  of  lines  (l-D's)  .  The  data  j,  p^ ,  q 
describe  the  lines  as  follows: 

j  -  the  index  of  j  •  1,...,N1, 

(p , ,q.)  m  the  index  numbers  of  the  endpoints  of  the  line 
g^  —  the  boundary  conditions  of  line  ,  where 
"0  ->  u^  free,  u2  free, 

1  ->  u^  fixed,  u2  free, 

2  ->  u^  free,  u2  fixed, 

8^-  -  3  «>  u^  fixed,  u2  fixed,  ► 

5  ->  u1  linear,  u2  free 

6  ->  u1  free,  u2  linear 

>.  7  ->  u ^  linear,  u2  linear.  ^ 


1  *,•  ff«  ’/ 


If  a  component  u^,  or  Uj  is  fixed  or  linear,  then  the  same  component 

must  be  fixed  at  the  endpoints  P°  ,  and  P°  .  Furthermore,  if  it  is 

qj  Pj 

fixed  then  the  two  values  given  at  these  endpoints  must  be  identical. 

pj  "  ±  1/Rj  S*-Sne<*  reciprocal  of  the  radius  of  the  arc 

segment  Pj  .  The  orientation  determines  the  sign  of  the  radius  as  shown 
in  Figure  1.4. 


“i  '  *1/R3 
Figure  1.4.  Orientation  of  the  arcs. 


N2  is  the  number  of  2-D  domains.  The  2-D  domains  are  prescribed  by 
the  4  corner  point  indices  if  N2  is  not  preceded  by  a  minus  sign,  and  by 
the  4  boundary  line  indices  if  N2  is  preceded  by  a  minus  sign. 

2 

The  first  input  in  each  line  after  N2  is  k  *  the  index  of  P^,  k=l,...,N2. 

12  3  4 

When  no  minus  sign  precedes  N2,  ir^,  n^,  ir^,  ■  the  index  numbers  of  the  4 

2 

cornerpoints  of  P^  ,  in  the  order  shown  in  Figure  1.5. 

c  4 


Figure  1.5.  Order  of  the  4  cbrnerpoint  indices  for  a  2-D  domain. 


The  user  should  be  careful  not  to  change  this  orientation.  For  example, 
the  Indices  cannot  be  input  in  the  order  as  displayed  in  Figure  1. 6. 

6 


„  *  J*  . 


•  *  _-  *_*  _«  V  .  , 


WRONG! 


Figure  1.6.  Illegal  ordering  of  the  4  cornerpoint  indices  for  a  2-D  domain 


12  3  4 

When  a  minus  sign  precedes  N2,  ir^,  ir^,  *  the  index  numbers  of  the  4 

2 

boundary  lines  of  V ^  in  the  order  shown  in  Figure  l.Z. 


Figure  1.7.  Order  of  the  4  boundary  line  Indices  for  a  2-D  domain. 
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Chapter  II.  Bilinear  and  Error  Matrices  Input 
2.1.  Introduction 

In  the  previous  section,  it  was  seen  that  FEARS  works  on  a  substructured 
domain,  which  allows  a  simple  input  and  a  more  accurate  representation  of  the 
geometry.  The  substructuring  has  another  advantage,  namely,  it  also  allows 
for  a  different  set  of  material  properties  (bilinear  matrices)  to  be  described 
on  each  2-D  domain. 


2.2.  Variational  Form 

The  variational  problem,  solvable  by  FEARS,  can  be  put  in  the  general 
form:  find  U  in  Mp  such  that 


£  uw\[i] ♦ v\[i]  -  [*r*  -  vclD}*  * 


N1  r  T  N2  f  r|* T  \  N1  f 

+  Z  vVu  ds-  Z  i  \4£\  D(x,y)  +  V1E(x»y)fdx  dy  +  Z  \  e .  (s)V  ds  , 
j-1  jjL  j  i*lJ  U3ZJ  J  j-lL  1 

vi  tr± 

for  all  V  in  ,  where  M^  and  M2  are  the  appropriate  finite 
dimensional  trial  and  test  spaces  (see  FEARS  Mathematical  Description) . 
Here,  we  have  used  the  notation 


U 


3U 

9Z 


r”3u^-i 

3x“ 
3u2 

dx~ 

ay 

3«2 
3y  J 


jV 

,  and  analogously  for  and  V  . 


Also, 

A^  is  a  A  x  A  symmetric  matrix, 
is  a  2  x  A  matrix, 

8 


.*• 

-rf-j 


.  fa.  A,  - 


C1  is  a  2  X  2 symmetric  matrix 


D^(x,y)  Is  a  4x1  vector  valued  function,  and 

^(Xjy)  Is  a  2  x  1  vector  valued  function,  for  1  _<  1  £  N2  . 

Yj  is  a  2  x  2  symmetric  matrix,  and 

Cj(s)  is  a  1  x  2  vector  valued  function  for  1  j  £  N1  . 


2.3  Error  Norms 


The  error  is  approximated  in  the  norm 


2P 


,  where  * 


3U 

az 

T 

r  - 
3U 

az 

fi 

L  dx  dy 

3-1  d2'* 
L  j 

1 

l/2p 


(2.1) 


The  value  p  is  a  global  value  which  is  input  before  we  input  the 
geometry  and  bilinear  forms.  The  values  P^  are  to  be  specified  for  each 
2-D  domain.  Normally  P^  *  Pj  "  “  PN2  *  P  »  which  will  correspond  to  a 
global  L2p  estimate.  The  corresponding  norm  (2.2),  is  used  when  p±  »  0 


(2.2) 


H  u* 


4x4  matrix,  and 


usually,  (Ag)  -  (A)j  . 


2.4  General  Format 

The  format  for  the  input  of  the  bilinear,  error,  and  output  matrices  is 
as  follows: 

*That  [TFTTTo  is  a  norm,  and  not  a  seminorm  for  the  error  e  ,  follows 
from  the  fact  cRat  if  |J|e[|L  •  0  ,  the  e  must  be  constant.  Since  the 
approximate  solution  is  bilinear,  the  exact  solution  must  be  bilinear  as 
well.  Because  of  quasi-optimality  of  the  finite  element  solution,  e*0  . 


NB 

Oj |  kj,  kj  , • 

...  k 

i  Bj  >  Yj »  Bj» 

eJ 

(A)j 

(if  Oj  i  0) 

(if  i  0) 

(o} 

(if  Yj  0) 

(if  fij  *  0) 

m3 

(if  tj  t  0) 

(**)i 

(% 

(S)j 

ij 


8I*  ei 


(Bilinear  and  Error 
Matrices) 

(This  packet  is  repeated 
NB  times ,  j«l, . . .  ,N1‘.) 


(Output  Matrix) 

(Line  Integration 
Matrices) 

(This  packet  is  repeated 
NL  times, i-1, ... ,NL.) 


2.5  Description  of  Bilinear  and  Error  Matrices 

NB  is  the  number  of  different  packages  of  bilinear  and  error  matrices. 
Each  different  package  must  then  be  listed.  In  each  package 
j  -  the  index  number  of  the  package, 

Oj  *  the  number  of  2-D  domains  for  which  the  package  of  matrices  applies 

k.,...,k  ■  the  indices  of  the  2-D  domains  for  which  the  package  applies. 

1  nj 

10 


°j*  Yj »  Cj »  are  integers  which  indicate  whether  or  not  the 

matrices  (A)^  ,(B)j  ,(C)j  ,(D)j  ,(E)^  are  aero,  in  which 
case  its  corresponding  input  line  is  not  present. 


0  ■>  (A) j  is  zero,  no  (A) ^  input  line, 

1  ->  (A) j  1*  0  ,  (A)^  line  is  present, 

0  «>  (B)^  -  0  ,  no  (B)^  input  line, 

1  ->  (B)j  j*  0  ,  (B)^  line  is  present. 


C  0  ->  (A)  j  is 
\l  ->  (A)j  1* 

CO  ->  (B^  - 

\l  ->  (B)j  j* 

f  0  “>  <c)j  - 

i  ll  ->  (C)  .  * 


0  ,  no  (C)j  input  line, 

0  ,  (D) ^  line  is  present. 


10  ■>  (D).  »  0  ,  no  (C).  input  line, 

1  ■>  (D)j  ^  0  ,  D  is  constant  and  is  defined  in  the  (D)^  line, 

-1  *>  (D)j  j*  0  ,  D  is  coordinate  dependent  and  must  be  defined 

in  a  subroutine;  dummy  values  must  be  supplied 
in  the  (D)^  line, 

(0  «>  (E)  ■  0  ,  no  (E).  input  line, 

1  ■>  (E)^  +  0  ,  E  is  constant  and  is  defined  in  the  (E)j  line, 

-1  ■>  (E)j  j*  0  ,  E  is  coordinate  dependent,  and  must  be  defined 

in  a  subroutine  and  some  dummy  values  must 
be  input  in  the  (E)^  line. 

Note:  If  ■  -1  or  ■  -1  a  subroutine  must  be  defined.  This  is 
described  in  Section  2.7. 

The  (A)j  line  should  contain  the  coefficient  values  of  the  A  x  A 
matrix  (A) ^  in  the  following  order, 


11 


•O.  O  O  . 


4  |  10 

~7  ~  13 


in  free  format.  The  appropriate  matrix  A  ,  when  using  plane  stress  or 

plane  strain  elasticity,  is  given  in  the  Appendix  (A. 2). 

T 

The  (B)j  line  should  contain  the  matrices  B  and  B  ,  where 


[V  *2] 


and  B 


bU>  b(J> 

bll  b12 


i  h(j)  K(j) 
b21  b22 


j  -  1,2  . 


The  matrices  should  be  input  in  the  form 


l  T 
Bi 


,  the  components 


of  this  matrix  being  input  in  the  same  order  as  with  (A) ^  . 

The  (C)^  line  should  contain  the  coefficient  values  of  the  2x2 
matrix  (C) ^  in  the  order 


C  .1 . 


The  (D)^  line  should  contain  the  coefficients  of  the  2x2'  matrix  (D) 
The  (E)  j  line  should  contain  the  entries  of  the  1x2  vector  (E)  ^  in 


the  order 


el*  *2 

The  (Ag)  line  should  contain  the  entries  for  the  4x4  error  matrix 
(Ag  )  in  the  same  order  as  the  entires  for  the  matrix  (A)^. 

The  (N  )  line  is  a  line  where  we  input  4  parameters  for  the  run  , 

C  j 


m 


Pj  ■  the  p  norm  for  the  domains  (see  2.1). 

If  %  <  p.  <  •  then  our  indicators  will  be  based  on  an  L_ 

i  -jpj 

estimate  over  these  domains.  The  value  Pj»0  ,  corresponds  to 
an  L9  based  error  estimate. 

Tj  -  the  weight  for  the  residual  part  of  the  error  indicator  which  is 
computed  through  an  integration  instead  of  a  derivative  jump. 

See  FEARS  Mathematical  description  for  more  details  on  this.  For 
most  cases  set  r*l  . 

Wj  *  a  parameter  which  was  formerly  used  in  the  residue  computation, 
^faysjjet^^t^jl  * 

Xj  ■  a  free  parameter.  Input  anything  here.  This  value  is  overwritten 
in  the  program. 


At  the  end  of  each  packet  we  must  input  the  matrix  for  output  generation 
(S) j .  (S)^  is  a  5  x  6  matrix  and  should  be  input  in  the  order 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

In  the  output  FEARS  will  compute  the  product 


run 

L  3Z 

LuJ 


for  2-D  domain  j  . 


The  appropriate  formulas  S^  ,  needed  to  yield  the  five  values  U^,  Uj* 

a  •  a  »  oOTi  for  the  plane  strain  and  plane  stress  assumptions  of  elasticity 
xx  yy  xy 

are  given  in  the  Appendix  (A. 2). 


2.6  Lin*  Integration  Boundary  Conditions 


NL  is  the  number  of  different  packages  of  line  integral  matrices.  These 
line  integrals  arise  from  boundary  conditions  in  which  traction  or  pressure 
forces  are  present.  These  forces  can  be  either  globally  defined ,  if  the  force 
is  in  a  fixed  global  direction,  or  locally  defined,  if  for  example,  the  force 
is  normal  to  the  boundary  as  with  a  hydrostatic  force.  The  vector  e  should 
give  the  magnitude  of  the  force  in  the  x  and  y  directions,  if  the  force 
is  global,  and  in  the  tangential  and  normal  directions,  if  the  force  is  locally 
defined.  The  2x2  matrix  y  is  present  whenever  the  force  depends  on  the 
displacement  (eg.  a  spring  force). 

Mathematically  the  boundary  conditions  allowed  are  of  the  form 


[4 


+  yU  ■  e  ,  where 


.30 

an 


.  rnx0"y°l  [jo 

0  nx  0  n  TZ 


,  where 


(n  ,n  )  is  the  outward  normal  to  the  boundary  and 
x  y 

2.2. 


3u"| 

3Zj 


is  defined  in  Section 


In  terms  of  stresses  the  boundary  conditions  are  of  the  form 
T  +  yU  ■  e  , 


where 


ft  M,  ♦  Ml 
L°xy  Mx  +  V  "yj 


-  [force  in  x  direction! 

-  [force  in  y  directionj 

When  the  forces  are  specified  locally  they  are  transformed  into  global  forces, 
which  change  in  direction  and  possibly  magnitude,  with  respect  to  the  arc 
length  of  the  boundary  line. 


I9WM 


After  NL  is  specified,  each  of  the  NL  packets  of  line  integrals  must 
be  specified.  In  the  line  following  NL,  we  must  input 
i  ■  the  index  number  of  the  package, 
n^  ”  the  number  of  1-D  domains  to  which  this  package  applies, 

*  the  indices  of  the  1-D  domains  to  which  the  package  applies. 

The  next  line  contains  indicators  g^,  e  where 

g^  *  ”0  ,  (y)±  ■  0  but  dummy  values  in  the  data  line  for  (y^ 

*  must  be  input, 

j.1  ,  (y)*  *  o  , 

ei  »  -1  ,  (c)A  y  0  and  the  values  of  (e^  describe  a  local  force 
on  the  boundary, 

0  ,  (e)  .  -  0  but  dummy  values  for  the  data  line  (e).  must 
«  1  1 

be  input, 

X  >  (c)j^  4  0  and  the  values  of  (e^  describe  a  global  force 
on  the  boundary. 


The  (Y)^  line  should  contain  the  entries  of  the  2x2  matrix  Y  in  the 


order 


(l  l) 


The  (e)^  line  should  contain  the  four  entries  e^(p),  e2(p),  e2(q),  e2(q) 
where  p  and  q  are  the  endpoints  of  the  boundary  line.  If  the  force  is 
described  globally  then  is  the  force  in  the  x  direction  and  c2  is 

the  force  in  the  y  direction.  If  the  force  is  described  locally  then 
and  e2  denote  the  forces  in  the  tangential  and  normal  directions 
respectively,  to  the  boundary. 

For  example,  suppose  we  wish  to  prescribe  a  force  of  magnitude  M  which 
makes  an  angle  of  a  with  the  tangent.  Then  the  appropriate  values  for  c 


*  *  *rn  *4 m  Vw  •  m.  fn  .»  v  ", 4  S  .  rVch  »"  «•“  «"  n  ,  •*  i  •"  .  »  ,  .  ,  •’  «  ,  •  J  •  ,  •  • 


should  be 


M  cos  a,  M  sin  a,  M  cos  a,  M  sin  a. 

Furthermore,  the  angle  and  magnitude  can  be  changed  linearly  (with  respect  to 
the  arc  length)  from  a  and  8  and  Mq  to  ,  respectively  by  Inputting 
Mq  cos  a,  Mq  sin  a,  cos  8,  sin  8  for  e  . 


2.7  Subroutines  for  the  functions  E(x,y),  D(x,y) 


If  or  ,  the  indicators  for  the  vectors  D  and  E  ,  are  equal 


to  -1,  we  must  supply  a  fortran  subroutine  to  define  these  functions. 

The  subroutines  will  read  in  the  values  (x(l),x(2)),  the  global  (x,y) 
coordinates  of  a  point  and  then  compute  and  return  the  values 
S(l)  “  the  exact  solution  u(l)  (if  known), 

S(2)  -  ^  the  exact  solution  u(2)  (if  known), 

D(l) 

I 

the  4  scalar  functions  comprising  the  components  of  the 


D(2) 

D(3) 

D(4) 


E(l) 


vector  D  , 


E(2) 


■ 


the  2  scalar  functions  comprising  the  components  of  the  vector  E, 


W(!)  -  j~D(l)  + 


the  derivatives  needed  for  the  residue  computation 


#0)1 

W(2)  -  j^>(2)  +  j 

The  subroutine  actually  has  4  entry  points,  returning  the  values  S,  D,  E, 
and  W,  respectively.  The  calling  sequence  is 


SUBROUTINE  ZPMTRU(X,S) 

DIMENSION  X(2),  S(2) 

(returns  the  exact  solution  S(l),  and  S(2), 


COMMON/ FUPAR/NSUP ,FUP(12) 

(Common  block  of  function  parameters.  NSUP-the  number  of  parameters, 
and  FUP(12)“the  function  parameters — as  dimensioned  NSUP  <_  12.  These 
parameters  are  input  at  the  onset  of  the  program  (see  Section  5.2.) 

ENTRY  DMATX(X,D) 

DIMENSION  D(4) 

(returns  the  components  D(l),  D(2) ,  D(3),  D(4))» 

ENTRY  EMATX(X,E) 

DIMENSION  E(2) 

(returns  the  values  E(l) ,  and  E(2)) 

and 

ENTRY  DMATXD (X  ,W) 

DIMENSION  W(2) 

(returns  the  values  W(l),  and  W(2))  . 

Tills  subroutine  must  be  compiled  and  mapped  with  the  main  program.  This 
is  briefly  discussed  in  Chapter  5. 


Chapter  III.  Commands  and  Strategies 
3.1  Introduction — Overview  of  Strategies 

Before  describing  in  detail  the  various  options  and  commands  available 
to  the  user,  we  give  a  brief  overview  of  how  FEARS  operates.  We  hope  this 
will  help  clarify  what  each  instruction  actually  does. 

Let  us  assume  that  the  geometry  and  bilinear  matrices  have  been  input. 

An  initial  subdivision  is  then  performed  by  the  program,  and  a  solution  is 
obtained  on  this  initial  mesh.  Error  indicators  are  then  computed  from  the 
solution  values.  A  status  message  or  REPORT  is  then  printed  out  for  this 
initial  solution. 

After  this  initial  step,  the  user  has  many  options  for  subsequent  sub¬ 
dividing  (refining)  and  resolving.  In  FEARS  this  iterative  process  of 
subdividing  and  resolving  is  continued  until  either  the  solution  obtained 
is  within  a  specified  tolerance  of  the  true  solution,  or  the  user  runs  out 
of  computer  resou'  :.es  (money,  time,  or  storage).  Each  REPORT  message  contains 
the  approximate  relative  error  as  part  of  its  output. 

Ideally  we  would  like  to  employ  some  optimal  strategy  which  will  obtain 

for  us  the  desired  accuracy  with  the  least  computer  expense.  On  -  one  hand , 

an  "optimal  mesh"  is  always  desired,  that  is,  a  mesh  which  will  yield  the 

smallest  error  in  the  solution  for  a  fixed  number  of  degrees  of  freedom.  On 

the  other  hand, we  would  not  like  to  spend  too  much  money  in  order  to  maintain 

an  optimal  mesh  at  each  level.  For  example,  it  would  be  very  expensive  and 

wasteful  if  we  subdivided  only  one  or  two  elements  at  a  time  and  then 

recomputed  the  entire  solution  after  each  subdivision.  Thus,  even  though  we 

may  get  a  better  mesh  by  subdividing  only  one  element  at  a  time,  it  would  be 

a  better  strategy  to  refine  a  larger  number  of  elements,  even  though  the 

mesh  obtained  may  be  only  "nearly  optimal".  It  has  been  shown  *(1)  that 

*(l)Babuska,  I.,  and  Rheinboldt,  W.  Analysis  of  finite-element  meshes  in  |R*. 
Math  Comput.  33  (1979),  435-463). 


(at  least  for  one  dimensional  problems)  a  mesh  which  deviates  slightly  from 
the  optimal  mesh  is  nearly  optimal  in  the  sense  that  the  solution  with  this 
mesh  is  nearly  as  accurate  as  the  solution  on  the  optimal  mesh. 

Furthermore,  one  may  ask  if  it  is  necessary  to  resolve  the  problem 
globally  after  each  subdivision.  For  example,  there  may  be  circumstances 
where  only  one  or  two  elements  will  get  subdivided,  or  all  the  elements  to 
be  subdivided  are  concentrated  in  one  region.  Perhaps  it  would  be  acceptable 
to  resolve  the  problem  locally — either  within  each  element,  or  only  within 
a  2-D  domain  where  some  subdivision  has  occurred.  The  FEARS  program  allows 
us  these  options. 

For  example,  we  can  specify  that  when  some  element  gets  subdivided,  all 
previously  obtained  solution  values  will  remain  fixed,  and  a  new  solution 
will  be  obtained  only  for  the  node  created  at  the  center  of  the  element  by 
the  subdivision.  This  is  referred  to  as  a  SHORT  path  solution.  Error 
Indicators  are  recomputed  for  only  the  4  new  elements.  SHORT  path  solutions 
are  fast  and  cheap  and  are  recommended  when  only  a  few  elements  are  to  be 
subdivided. 

Likewise  we  may  specify  some  set  of  2-D  domains,  for  example,  only 
those  where  subdivision  has  occurred,  and  then  obtain  new  solution  values 
only  for  those  2-D  domains.  The  boundary  conditions  on  the  boundary  of 
the  subdomain  are  taken  to  be  fixed,  with  displacement  values  determined 
from  the  previous  solution.  This  type  of  solution  path  is  in  between  a 
SHORT  path  and  full  solution  in  both  expense  and  accuracy. 

Although  the  user  can  control  the  refinement  procedure  by  specifying 
which  elements  are  to  be  subdivided  before  each  solution  path  FEARS  also  has 
a  built  in  recommended  refinement  strategy.  This  strategy  is  enacted  through 
the  command  AUTO  (short  for  automatic) .  Each  time  this  command  is  given, 


all  those  elements  with  error  indicators  larger  than  some  computed  threshold 
value  (see  Appendix  A3)  are  subdivided  and  a  new  solution  is  computed.  The 
type  of  solution  path  performed  must  be  supplied  by  the  user.  For  example, 
(AUTO/1)  will  refine  and  recompute  the  solution  globally,  and  (AUTO/4)  will 
refine  and  then  perform  a  SHORT  path  solution  on  each  refined  element  only. 

For  many  problems,  particularly  those  in  which  there  are  no  singularities, 
a  sequence  of  (AUTO/1)  commands,  is  the  best  strategy  for  obtaining  an  accurate 
solution  cheaply.  However,  when  solving  problems  with  singularities,  it  is 
often  the  case  that  an  AUTO  command  will  refine  only  one  or  two  elements. 

If  the  mesh  already  has  a  large  number  of  elements,  producing  a  new  global 
solution  in  this  case  is  not  only  costly,  but  also  unnecessary.  In  this 
case  we  would  prefer  performing  a  SHORT  path  solution  with  an  (AUTO/4) 
instead  of  a  new  global  solution  with  (AUTO/1).  If  the  program  is  being 
run  interactively,  then  the  user  can  decide  which  type  of  solution  should 
be  performed,  since  after  each  REPORT  the  approximate  number  of  elements 
to  be  subdivided  by  the  next  AUTO  command  is  printed. 

Unfortunately,  if  the  program  is  being  run  as  a  batch  job  there  is  no 
a-priorl  way  to  determine  when  an  (AUTO/1)  or  (AUTO/4)  should  be  performed. 
FEARS  also  has  the  ability  to  make  this  choice  automatically  with  the  ITER 
(ITERATE)  command.  The  ITER  command  performs  n  solution  paths  composed 
of  (AUTO/1)  or  (AUTO/4)  commands,  the  choice  depending  on  whether  the  number 
of  elements  in  the  new  mesh  is  a  certain  percentage  over  the  number  in  the 
old  mesh  (after  the  last  (AUTO/1)).  This  cut-off  percentage  must  be  supplied 
by  the  user,  but  a  30%  increase  is  recommended. 

The  user  also  has  control  over  what  information  is  printed  out,  after 
a  solution  path  is  performed.  For  example,  with  the  PRINT  command,  informa¬ 
tion  about  the  solution  at  the  nodes  and  a  list  of  elements  with  their  a- 
posteriorl  error  indicators  may  be  printed.  The  OUTPUT  command  wi31  give  a 
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list  of  stresses  and  solution  values  at  the  center  of  each  element. 

Two  other  useful  commands  are  the  DUMP  command  and  the  CHANGE  command. 

The  DUMP  command  will  cause  all  information  about  the  problem,  eg.  geometry, 
bilinear  matrices,  solution  values,  and  data  structures  for  the  mesh,  to 
be  saved  on  a 'permanent  file.  The  user  may  then  restart  the  program  where 
he/she  left  off  at  any  future  time. 

The  CHANGE  command  causes  small  changes  in  the  initial  problem.  This 
is  useful  if  one  is  interested  in  the  effect  of  perturbing  either  the 
geometry  or  material  properties  of  the  prolem.  In  this  case,  the  refined 
mesh  for  the  original  problem  will  be  almost  topologically  equivalent  to 
the  refined  mesh  for  the  next  problem.  Therefore,  instead  of  using  a  lot 
of  computer  time  by  iteratively  subdividing  and  resolving  for  the  new  problem, 
the  refine4  mesh  for  the  original  problem  could  be  used  and  a  final  solution 
obtained  immediately. 

Now  that  the  general  format  has  been  presented,  we  describe  the  user 
commands  in  detail.  The  computer  will  always  acknowledge  that  it  is  ready 
to  receive  a  user  command  by  printing  the  line 

****  COMMAND: 

3.2  PRINT  Command 

The  PRINT  command  is  designed  to  print  out  information  about  the 

points  -  D ®  ,  j«l,...,NO,  the 

lines  -  ,  j«l,...,Nl,  and  the 

2 

2-D  domains  -  ,  j«l,...,N2  . 

The  format  for  the  PRINT  command  is 

PRINT 

a,b,c 
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K.  ^ 


h 


The  value,  a,  determines  the  dimension  of  the  domains  to  be  printed. 
0  ■>  print  about  Z^’s 

1  *>  print  about  d\* s 


2  ->  print  about  ZJ*' s 

0  1  2 

-1  *>  print  about  D.  's,  D, *s  and  0. 's  . 
W  3  J  J 


The  value,  b,  determines  which  index  k  of  Z^Ca^-l)  is  to  be  printed. 


{>  i  -  k , 
-1  pri 


k  ,  print  about  Z?*  only, 
print  all  Z^'s  . 


The  value,  c,  determines  how  much  information  is  to  be  printed. 

* 

0  print  only  headings 

1  only  print  data  about  the  nodal  points  of  the  2-D  domain(s) , 

c  ”  « 

2  only  print  data  about  the  elements  of  the  2-D  domain (s). 

3  print  all  information 


For  example,  the  command 
PRINT 
-1.-1.3 

will  cause  all  information  about  each  0-D,  1-D  and  2-D  domain  to  be 
printed.  The  format  of  the  printed  data  and  its  interpretation  are 
discussed  in  detail  in  Chapter  4. 

After  the  execution  of  a  PRINT  command  the  computer  will  return  to 
the  command  mode  and  print  ****  COMMAND  . 


v.v.v.v;,'.' 
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[>  *.% 

K 


3.3  REPORT  Command 

The  REPORT  command  accomplishes  two  things.  First  of  all,  it  computes 
statistical  information  about  the  error  indicators  and  sorts  the  elements 
in  the  order  of  decreasing  indicators.  Thus,  a  REPORT  should  be  performed 
before  each  PRINT  command  to  ensure  that  the  elements  are  listed  in  order. 
Secondly,  it  prints  out  a  status  report  on  the  full  domain  giving  informa¬ 
tion  on  the  number  of  elements,  total  energy,  error  estimator,  percentage 
of  error,  number  of  elements  recommended  for  subdivision,  etc.  This  message 
is  described  in  detail  in  Chapter  4-Output. 

3.4  SUBDIVIDE  Command 

The  SUBDIVIDE  command  is  used  to  subdivide  elements.  It  has  the 
following  format: 

computer  ****  COMMAND: 


user 


SUBDIVIDE 


computer  SUBDIVISION,  2-D  PROCESS  INDEX: 


user 


computer  Q  -ALL ,  GT .  0 .  -CUT  OFF  VALUE , 

LT.O.  -  GIVE  ELEMENT  LIST 


user 


computer  f  ELT  LVL  R  LOCAL  COORD  ERR.  IND  PREV  ERRIND 


user 


ELEMENT  TO  BE  SUBDIVIDED: 


IF  %1  <  0 


computer 


ELEMENT  TO  BE  SUBDIVIDED: 
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user 


computer 


user 


computer 


usbr 


computer 


ELEMENT  TO  BE  SUBDIVIDED: 


.  0  (zero) 

SUBDIVISION  2-D  PROCESS  INDEX: 


0  »  ALL,  GT.O.  -  CUT  OFF  VALUE, 
LT.O.  -  GIVE  ELEMENT  LIST. 


IF  Xj  <  0 


user 


computer  SUBDIVISION,  2-D  PROCESS  INDEX: 


user 


(zero) 


computer  ****  COMMAND: 

The  values  give  the  indices  of  the  2-D  domains  in 

the  subdivision  will  occur. 


which 


The  values  X^,  X^,...,  determine  which  elements  get  subdivided. 

X^  ■  0  *>  subdivide  all  elements  in  2-D  domain  with 

<  an  error  indicator  >_  X^. 

* 

0  =>  give  an  element  list  and  ask  for  an  index  of  an 
element  to  be  subdivided.  The  program  will 
continue  to  ask  for  an  element  until  the  user 


returns  0  (zero). 

The  user  may  get  out  of  the  SUBDIVIDE  loop  with  a  0  (zero) . 


3.5  DOMAIN  Command 

The  DOMAIN  command  is  used  to  set  up  the  subdomain  on  which  a  solution 
path  is  to  be  performed.  The  subdomain  defined  may  be  any  subset  of  the 
D*'b  . 
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The  format  Is  either 

(I)  DOMAIN 

0 

or 

(II)  DOMAIN 
N 

k.  * • • • ik  • 

1  N 

In  case  (1)  the  subdomain  defined  is  the  full  domain,  and  in  case  (ii)  the 
subdomain  consists  of  the  N  2-D  domains  •  After  the  command 

DOMAIN  is  returned  to  the  computer,  the  message 
SUBSET  OF  2-D'S:  NUMBER  OF  2-D'S  (0  -  ALL) 
vill  be  printed  by  the  program. 

3.6  LONG  Command 

The  LONG  command  will  obtain  new  solution  values  for  each  node  in  the 
subdomain  specified  in  the  DOMAIN  command.  If  the  subdomain  defined  with 
DOMAIN  is  a  proper  subdomain  of  the  full  domain  then  any  points  (D^'s) 
which  are  external  to  this  subdomain,  but  are  Internal  in  the  full  domain 
are  considered  as  fixed  with  the  values  prescribed  by  a  previous  solution 
path.  Also,  error  indicators  of  elements  which  are  not  in  the  subdomain 
are  not  recalculated. 

The  format  for  this  command  is  simply 

LONG  . 

3.7  SHORT  Command 

The  SHORT  command  performs  a  subdivision  and  a  short  path  solution 
on  each  element  specified.  When  an  element  is  specified,  it  is  sub¬ 

divided  into  4  sub-elements,  and  solution  values  are  obtained  for  only 
the  new  node  formed  in  the  middle  of  e^  .  This  solution  is  computed  by 
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solving  the  problem  on  the  4  sub-elements  with  linear  boundary  conditions 
on  the  boundary  of  e^  determined  by  the  4  cornerpoint  solution  values 
previously  determined  on  e.  . 


Figure  3.1.  Short  path  subdivision  of  element  e^. 

In  the  figure  we  solve  for  the  circled  node  in  the  center  having  the 
previously  obtained  solution  values  at  1,  2,  3,  and  4  and  the  values  at 
5,  6 ,  7,  and  8  are  obtained  thrpugh  interpolation.  The  new  error  indicators 
for  |^he  four  new  elements  created  are  computed  in  the  following  way. 

Let  F  denote  the  father  element  (to  be  subdivided)  and  S^,  Sj,  S^, 
and  the  four  sons  generated  from  subdividing  F.  Let  E(F)  denote  the 
error  indicator  for  the  father  element,  and  P(F)  the  predicted  error 
indicator  for  the  four  sons  of  the  father  (see  Figure  3.2).  Appendix  A. 3 
describes  how  this  predicted  error  indicator  is  computed. 


Figure  3.2.  A  father  element  subdivided  into  four  sons. 


We  then  compute  the  error  for  each  son  E^)  by  the  formula 
E(Sj)  -  min(aE(F)  ,  P(F)),  where 

a  m  .55  by  default,  but  can  be  changed  by  the  user  with  the  SHEF  command 


\ *v.  *.*  *.  * »»%.  *»* . *j5‘. 


(see  Section  3.16) (  or  at  the  start  of  the  prograa  (see  Chapter  5). 

The  format  for  the  SHORT  command  is  similar  to  the  SUBDIVIDE  command 


s 

computer 

****  COMMAND 

user 

SHORT 

computer 

2-D  DOMAIN  PROCESS  INDEX: 

>>• 

*v% 

user 

J1 

computer 

ELT  LVL  R  LOCAL  COORD  ERR 

el  A1  rl 


yl  *1 


computer  ELEMENT  TO  BE  SUBDIVIDED  AND  SOLVED 


computer 

user 

computer 

user 


ELEMENT  TO  BE  SUBDIVIDED  AND  SOLVED 
@EOF 

2-D  DOMAIN  PROCESS  INDEX: 


user 

computer 

user 

computer 


@EOF 

2-D  DOMAIN  PROCESS  INDEX: 
@EOF 

***  COMMAND 


3.8  AUTO  Command 

The  AUTO  command,  short  for  automatic,  is  a  powerful  and  useful  user 
command.  It  performs  the  subdivision,  sets  up  the  subdomain,  performs  a 
solution  path,  and  prints  a  new  report. 
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The  format  for  the  AUTO  command  is 
AUTO 
J 

X  (if  J  <  0) 

If  J  >  0  each  element  having  an  error  indicator  which  is  greater 
then  or  equal  to  a  computed  threshold  value  is  subdivided. 

If  J  <  0  another  value  X  must  be  input  and  all  elements  having  an 
error  indicator  greater  than  or  equal  to  X  will  be  subdivided. 

The  value  of  J  determines  which  type  of  solution  path  to  take. 

J  ■  ±  1  ■>  Obtain  a  new  solution  for  the  full  domain. 

J  ■  ±  2  m>  Obtain  a  new  solution  for  the  subdomain  composed  of  those 
2-D  domains  where  subdivision  occurred. 

J  ■  ±  3  ■>  Obtain  a  new  solution  for  each  2-D  domain  individually 
where  subdivision  occurred. 

J  «  i  4  »>  Perform  a  short  path  solution  for  each  element  subdivided. 
Thus,  if  we  wish  to  subdivide  all  elements  of  our  full  domain  and 
resolve  the  problem  globally  we  could  use  the  command 
AUTO 
-1 
0. 

The  same  effect  could  be  accomplished  with  the  sequence  of  commands: 
SUBDIVIDE 


DOMAIN 


0 

LONG 

REPORT 

When  using  the  AUTO  command,  it  is  important  to  understand  how  the 
threshold  value  is  computed.  Appendix  A. 3  contains  a  detailed  explanation. 
Actually  we  compute  a  threshold  for  each  2-D  domain  and  then  use  a  global 
threshold  being  the  maximum  of  all  the  2-D  domain  thresholds.  All  elements 
above  the  global  threshold  are  then  subdivided  in  the  AUTO  command  (if 
J  >  0). 

3.9  ERROR  Command 

The  ERROR  command  recalculates  all  error  indicators  using  the  solutions 
last  obtained  and  generates  a  new  report.  It  should  be  used  after  repeated 
short  path  solutions  in  order  to  obtain  more  accurate  error  Indicators. 

3.10  DEBUG  Command 

The  DEBUG  command  was  initially  used  as  a  debugging  aid,  and  this 
option  is  still  available.  However,  the  user  may  find  it  more  important 
in  its  ability  to  print  a  list  of  subdivided  elements  and/or  its  control 
of  the  output. 

The  format  for  the  DEBUG  command  is 
computer  ****  COMMAND: 

ueer  DEBUG 

user  IPR(l)  ,  IPR(2) , . . . ,IPR(8) 

Where  IPR(l)  , . . .  ,IPR(8)  deal  with  the  following: 


\  s 


N 


n 

* 

:•! 


VS 

•.n 


R 

!►  ■' 


r  » 


f? 


C 


IPR(l)  -  subdivision  element  list,  short  path  debug* 

IPR(2)  -  debug  error  calculation, 

IPR(3)  -  debug  matrix  assembly* 

IPR(4)  -  debug  matrix  decomposition*. 

IPR<5)  -  debug  matrix  solution, 

I?R(6)  -  debug  back  substitution* 

IPR(7)  -  echo  input, 

1PR(8)  -  automatic  elemental  output  control. 

The  values  IPR(l) . . .  IPR(8)  should  be  specified  as  follows: 


IPR(l) 


r  0 


1 


do  nothing, 

print  a  list  of  the  elements  which  were  subdivided 


2 

-k 


IPR(J)  - 


IPR(7)  - 


0 
1 
12 

{: 


IPR(8)  - 


before  each  solution  path* 

print  out  short  path  solution  debugging  Information 
record  the  subdivision  element  list  on  file  with  FORTRAN 
unit  number  k  .  File  k  must  have  been  properly 
assigned  to  the  run. 


do  not  print, 
print  only  summary* 
print  all  information, J 
do  not  echo  input, 
echo  input, 

no  print*  no  file  write, 
print*  no  file  write* 


(For  j  -  2,*.  ..,6.) 


2  no  print,  file  write* 

3  print,  file  write. 

The  parameter  IPR(8)  causes  the  listing  of  the  OUTPUT  information 
(see  Section  4.6)  to  be,  either  printed,  or  written  on  FORTRAN  unit  number 
17,  after  each  AUTO  command. 


The  program  also  asks  for  the  values  IPR(l) , . . . ,IPR(8)  at  the  beginning 
of  each  run,  and  the  DEBUG  command  offers  a  way  to  change  the  Initial  values 
of  these  parameters. 

3.11  OUTPUT  Command 

The  OUTPUT  command  will  have  the  effect  of  temporarily  changing  the 
parameter  XPR(8)  and  immediately  performing  the  file  write  and/or  print  as 
desired.  The  format  is 

computer  ****  COMMAND: 
user  OUTPUT 

user  n  (n-1,  2  or  3) 

This  effect  is  only  temporary.  If  a  new  AUTO  command  is  performed  the 
print  and/or  file  write  will  be  done  as  specified  by  either  the  last  DEBUG 
command,  or  by  the  initital  value  given  IPR(8)  at  the  start  of  the  program. 

3.12  DUMP  Command 

The  DUMP  command  allows  us  to  save  the  present  data  structures,  mesh, 
solution  values,  etc.  in  some  file  so  that  the  problem  can  be  restarted 
from  where  we  left  off  at  some  later  time.  The  format  is 
DUMP 
F 

where  F  is  an  integer  indicating  some  fortran  unit  number.  Since  units 
10-17  are  already  used  in  the  program,  units  18  and  up  may  be  used  here. 
These  files  should  all  be  initially  assigned  before  the  program  run,  and 
this  is  described  in  Chapter  5. 

3.13  RESET  Command 

The  RESET  command  will  restart  the  problem  from  the  time  just  before 
the  last  DUMP  was  executed.  The  format  is 


RESET 


F 

where  F  is  the  same  unit  number  used  with  the  DUMP  command.  If  you  were 
presently  running  another  problem  the  RESET  command  will  destroy  your  present 
run  unless  you  DUMP  it  onto  some  other  unit  number. 

3.14  ITERATE  Command 

The  ITERATE  command  causes  the  program  to  take  iteration  steps, 
composed  of  AUTO/J  command,  automatically,  with  built  in  termination.  These 
steps  will  be  either  LONG  solutions  (J-l)  or  SHORT  solutions  (J=*4) . 

The  format  for  the  ITERATE  command  is: 
computer  ****  COMMAND: 

user  ITERATE 

user  n,  8,  t,  a 

where  n(>0)  is  integer  valued,  8,  t,  a  (>0)  are  real  valued  inputs.  The 
value  n  determines  the  maximum  number  of  AUTO/J  commands  to  be  performed 
in  the  sequence. 

The  value  of  8  is  used  in  the  decision  strategy  to  determine  the 
value  of  J(«l  or  4).  The  decision  value  works  in  the  following  way.  Let 
EL  be  the  number  of  elements  after  the  last  LONG  solution  (AUTO/1).  Let  ES 
be  the  estimated  number  of  elements  in  the  mesh  when  the  next  refinement 
occurs.  Then  if  ES  <  (l+S)  EL,  the  AUTO/4  command  (SHORT)  is  performed, 
without  changing  the  value  of  EL.  Otherwise,  if  the  estimated  increase  in 
elements  is  8%  or  more,  an  AUTO/1  (LONG)  is  performed  which  also  changes 
EL  to  the  new  number  of  elements.  Exception  to  this  rule  is  at  the  last 
step  of  the  iteration  steps,  when  the  program  always  generates  an  AUTO/1 
(LONG  solution) . 

The  value  t  is  the  maximum  allowed  time  in  seconds.  The  program 
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assumes  that  the  time  needed  for  an  AUTO/1  (LONG)  step  is  linearly  dependent 
on  the  number  of  elements  E  ; 

T1  -  C  *  E  , 

where  the  factor  C  is  derived  from  a  previous  AUTO/1  step.  If  T1  and  the 
previously  accumulated  time  of  the  iteration  process  is  equal  or  greater  than  • 
the  given  t  value,  then  a  last  step,  AUTO/1,  is  performed  to  assure  a  LONG 
solution  before  returning  to  the  user's  next  command. 

The  value  a  is  the  required  relative  accuracy.  Again,  if  a  has 
been  achieved  in  the  sequence  of  AUTO/J  steps,  the  program  assures  that  the 
last  step  has  been  a  LONG  solution  before  returning  to  the  user's  next 
command. 

During  the  sequence  of  AUTO/J  steps,  it  may  happen  that  the  available 
data  storage  area  is  exhausted  due  to  the  increasing  number  of  elements.  The 
required  storage  areas  are  also  estimated  by  assuming  linear  dependency  on 
the  number  elements,  although  this  may  not  be  very  accurate. 

Thus,  an  ITERATE  command  causes  a  sequence  of  AUTO/J  steps,  where  the 
number  of  steps  is  determined  by  one  of  the  four  factors: 

(i)  n  ■  given  maximum  number  of  steps  (no  message  is  printed) 

(ii)  time  allowance  (t)  is  exhausted 

(iii)  accuracy  (a)  achieved 

(iv)  storage  is  exhausted 

If  one  of  the  last  3  cases  occurs,  the  appropriate  message  is  printed. 

3.15  CHANGE  Command 

The  CHANGE  command  allows  the  user  to  make  modifications  in  either  the 
bilinear  form,  geometry,  or  function  parameters,  while  keeping  the  mesh  generated 
from  solving  the  original  problem.  This  command  is  useful,  if  for  example,  the 
problem  of  optimal  design  is  of  interest. 


A  problem  is  solved  with  some  initial  geometry  using  the  FEARS  program 
and  an  optimal  mesh  is  created  for  this  problem.  It  may  be  of  interest,  for 
example,  to  determine  how  the  maximum  stress  is  altered  by  making  a  perturbance 
of  the  geometry.  With  the  CHANGE  command,  this  perturbed  problem  can  be  solved 
with  just  one  solution  pass,  using  the  final  mesh  of  the  original  problem. 

Suppose  that  the  final  mesh  of  the  original  problem  was  saved  on  file 
k  by  using  the  DUMP  command.  After  calling  the  FEARS  program  (see  Chapter  5) 
you  will  first  be  asked  to  supply  the  8  DEBUG  parameters  IPR(1)-IPR(8) . 

After  supplying  these  values  the  computer  will  respond: 
computer  PROGRAM  INPUT  DATA  IS  ON 

FILE  NUMBER  = 
to  which  you  should  input 
user  k 

k  is  the  file  number  on  which  the  old  mesh  was  stored. 

You  will  then  be  in  the  command  mode  as  the  computer  will 
respond 

****  COMMAND:  . 

To  use  the  CHANGE  command  the  user's  response  is 
user  CHANGE 

The  program  will  then  shift  back  into  the  input  mode  and 
ask  for  the  function  parameter^  ID  Number  of  the  problem, 
P-NORM  for  the  full  domain,  and  then  the  geometry  and 
bilinear  forms.  Here  the  data  should  be  prepared  with 
the  appropriate  changes,  and  input.  The  program  will  then 
immediately  obtain  a  solution  for  the  new  problem  using  the 
previously  saved  mesh.  A  REPORT  will  be  given  and  the 
program  will  return  to  the  command  mode. 
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3.16  SHEF  Command 


SHEF  stands  for  short  £ath  error  _factor.  This  command  allows  the  user 
to  prescribe  the  value  of  a  which  is  used  in  the  short  path  error  calcula¬ 
tions  (see  Section  3.7).  The  format  for  the  SHEF  command  is 
computer  ***  COMMAND 
user  SHEF 

user  a 

where  a  >  0  is  real.  If  this  command  is  not  used  the  default  value  for 
a  is  0.55. 

3.17  STOP  Command 

The  STOP  command  causes  the  termination  of  the  program  execution,  thus 
it  should  be  the  last  given  command: 
computer  ****  COMMAND 

user  STOP 

3.18  ERIT  Command 

The  ERIT  command  changes  the  way  one  of  the  error  terms  is  calculated 
for  elements  adjacent  to  internal  1-D  domains.  This  error  term  is  based 
on  the  difference  of  derivatives  of  the  solution  across  the  interfacing 
1-D  domain.  Initially  in  the  normal  (ON)  position,  the  difference  of 
derivatives  is  computed  form  the  two  adjacent  2-D  domains;  while  in  the  OFF 
position,  the  difference  of  derivatives  are  estimated  internally  in  the 
appropriate  2-D  domain.  The  format  of  the  ERIT  Command  is 
computer  ****  COMMAND 

user  ERIT 

which  causes  the  change  of  position,  from  ON  to  OFF,  or  OFF  to  ON.  This 
command  is  particularly  useful  when  the  problem  contains  an  interface 
separating  two  2-D  domains  with  vastly  different  material  properties. 
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3.19  MESH  Command 


The  MESH  command  Is  similar  in  purpose  to  the  CHANGE  command,  l.e.t  to 
get  solutions  for  a  slightly  changed  (In  geometry  or  material  constants) 
problem.  The  use  of  the  MESH  command  requires  that  the  subdivision  element 
list  of  the  original  problem  has  been  saved  by  the  DEBUG  command 
DEBUG 
— k | ... 

during  the  iteration  process  on  file  k  ,  where  it  is  a  Fortran  unit 
number  properly  assigned  to  the  run.  If  this  has  been  done,  one  can  use  the 
MESH  command  after  the  initial  input  (geometry,  bilinear  forms)  to  recreate 
the  saved  mesh  structure  for  the  modified  problem  by 
computer  ****  COMMAND 

user  MESH 

user  k 

Note,  that  the  MESH  command  only  recreates  the  mesh  according  to  the  data 
stored  on  file  k  and  does  not  obtain  a  solution  on  it.  Then  it  should 
be  followed  by  a  LONG  command  if  the  user  wishes  to  obtain  solution  on  the 
recreated  mesh. 

3.20  INIT  Command 

The  INIT  command  reinitializes  the  program  and  expects  initial  input 
as  a  new  start  without  asking  for  the  output  option. 
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CHAPTER  IV.  THE  OUTPUT 


4.1  Introduction-Global  and  Local  Coordinates 

Recall  that  FEARS  takes  a  unit  square  and  naps  it  onto  each  2-D  domain. 
All  computations  and  mesh  refinements  are  actually  computed  on  this  unit 
square  under  the  appropriate  transformation.  The  coordinates  (£*n)  on  the 
unit  square  are  called  the  local  coordinates  and  the  corresponding  values 
(Xty)"(x(£>n) ,  y(£»n3>  are  the  global  coordinates . 


Figure  4.1.  Local  and  Global  coordinates  for  a  2-D  domain. 
Inside  this  unit  square  all  refinements,  element  ordering,  etc.,  takes  place. 

4.2  Numbering  the  Mesh  of  a  2-D  Domain 

In  order  to  help  read  the  printout,  a  description  of  how  the  elements  and 
nodes  are  indexed  is  given.  This  will  be  a  great  aid  in  reconstructing  the 
mesh  and  hence  in  locating  elements  and  their  neighbors. 

The  initial  subdivision  is  numbered  in  the  following  way: 


Figure  4.2.  Numbering  for  initial  subdivision. 


The  nodes  and  elements  are  numbered  together— the  index  1  corresponding  to 
the  node  in  the  center  and  the  indices  2,  3,  4,  and  5  number  the  4  elements. 

Suppose  that  element  2  gets  subdivided.  Then  index  2  will  correspond 
to  the  node  formed  at  the  center  and  the  next  available  Indices  6,  7,  8,  and 

• 

9  will  be  used  to  number  the  4  new  elements  created  (see  Figure  4.3a).  Note 
that  the  two  circled  points  are  not  numbered.  These  points  are  irregular 
points  without  degrees  of  freedom,  whose  solution  values  are  obtained  through 
interpolation.  Generally,  regular  points  (nodes)  are  those  which  lie  at  the 
corner  of  4  elements  and  are  always  numbered.  Irregular  points  lie  at  the 
corner  of  only  2  elements  and  on  the  side  of  some  other  element  and  are  not 
numbered. 

Suppose  next  we  subdivide  element  3.  Then  index  3  is  a  point  and  the 
next  4  indices  10,  11,  12,  13  number  the  4  new  elements  created  (see  Figure 
4.3b).  However,  note  that  the  point  with  local  coordinates  (.25,  .5)  is  now 
a  regular  point  and  hence  gets  the  next  available  index  14. 

Finally  we  present  the  numbering  after  elements  4  and  5  get  subdivided 


in  Figure  4.3c 


A  listing  of  the  elements  subdivided  at  each  level  will  be  output  if 
the  parameter  IPR(l)  is  set  to  1  with  the  DEBUG  command  (see  Section  3.10). 


4.3  The  PRINT  output 


The  PRINT  command  will  cause  data  to  be  printed  about  the  points  (D^'s), 


1  2 
lines  (D^'s),  and  2-D  domains  (5j's)  . 


j 


4.3.1.  Printing  the  Points  (0-D  domains) 

The  format  for  the  data  about  the  0-D  domains  of  the  geometry  is  as 
follows: 


«<  0-D 

COORDINATES 

SOLUTION 

ERROR 

BDRY 

EXT 

1 

Xlyl 

U1  V1 

6U1  eVl 

bi 

ex. 

1 

2 

• 

x2  y2 
• 

U2  v2 
e 

eu2  ev2 

• 

b2 

• 

ex2 

• 

i. 

e 

• 

^O^NO 

e 

e 

^NO^NO 

e 

e 

eUN0eVN0 

• 

e 

bNO 

• 

eXNO 

The  information  under  the  heading 

0-D  -  gives  the  index  number  of  the  point, 

COORDINATES  -  gives  the  global  (x,y)  coordinates  of  the  point  as  specified 

by  the  geometry  input, 

SOLUTION  -  gives  the  computed  solution  values  of  the  point, 

ERROR  -  gives  the  error  between  the  computed  and  exact  solution  if 

the  exact  solution  is  known  and  supplied  in  the  subroutine 


ZMPTRU, 

BDRY  -  gives  the  boundary  conditions  of  the  point  as  specified  by 

the  geometry  input,  and 

EXT  -  indicates  if  the  point  is  internal  (0)  or  external  (1)  to 

the  full  domain. 
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4.3.2  Printing  the  lines  (1-D  domains) 


The  format  of  the  Information  about  the  line  with  index  J  is 

«<  1-D  INDEX:  J  B:  b,  E:  ex,  1/R:  p  , 

J 


FROM  aj  TO  bj  PTS:  p,  R-PTS:  rp 


j 

j 

P 

POINTS: 

PT 

R 

LOCAL 

GLOBAL 

COOR 

SOLUTION 

ERROR 

1 

• 

rl 

P1 

X1 

*1 

U1  V1 

8U1  CV1 

• 

P 

rP 

PP 

*P 

yp 

“P  VP 

eiip  eVp 

On  the  top  line  the  value  after 

B:  gives  the  boundary  condition  of  the  line  as  specified  by 

the  geometry  input, 

E:  indicates  if  the  line  is  internal  (0)  or  external  (1) 

1/R:'  gives  the  signed  reciprocal  of  the  radius  as  specified  by  the 

geometry  input, 

FR0M_T0_:  gives  the  two  indices  of  the  0-D  domains  forming  the  endpoints 
of  the  line, 

PTS:  gives  the  number  of  points  on  the  interior  of  the  line 

arising  from  the  mesh,  and 
R-PTS:  gives  the  number  of  regular  points. 

Next,  a  list  of  data  about  the  points  lying  in  the  interior  of  the  line 
arising  from  the  mesh  is  printed.  The  data  under 

PT  -  gives  an  index  of  the  point  in  the  order  it  was  formed, 

R  -  indicates  if  the  point  is  regular  (3)  or  not  (0,  1  or  2) , 

LOCAL  -  gives  the  local  coordinate  of  the  point  on  the  line — the 

0-D  aj  has  a  local  coordinate  0  and  bj  has  a  local  coordinate  1 


GLOBAL  COORDINATE  -  gives  the  global  (x,y)  coordinate  of  the  point, 
SOLUTION  -  gives  the  computed  solution  values  at  the  point,  and 
ERROR  -  gives  the  error  between  the  computed  and  exact  solutions  if 
the  exact  solution  is  known. 


4.3.3  The  2-D  Domain  Printout 

^The  format  for  the  printout  of  2-D  domain  J  is 


<«2-D  INDEX: J  CORNERS :  NUMBER  OF  POINTS :NP  ELEMENTS :  NE  ENERGY: 

ELT.SIZES:S1S2S3S4S5S6S7Sg,  SMALLER  SIZES : Sj, MAX. ADJ. RATIO: _ 

ERROR  INDICATORS,  TOTAL:  .MAX:  ,  MEAN:  ,  THRESHOLD :_ _ 


EAD-^  DISTR.  BELOW  MEAN  (_) : 
ING 


DISTR.  ABOVE  MEAN  (__):  _ _ _ _ 

DISTR.  ABOVE  PRED  (_):  _ _ _ _ 

TIME,  ASMS  DEC: _ ,  BCK:  THR: _ ,  TOTAL: 

STORAGE,  MEMORY: _ ,  AUXILIARY: 
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gives  the  4  coraerpoint  indices  as  specified  by  the 
geometry  input, 


NUMBER  OP  POINTS:  gives  the  number  of  regular  points  (NP)  of  the  mesh  in  the 


2-D  domain  J, 

ELEMENTS:  gives  the  number  of  elements  (NE)  in  the  mesh,  and 

ENERGY:  gives  the  energy  of  2-D  domain  J. 

In  the  next  line 

ELT.  SIZES:  gives  the  number  of  elements  in  the  domain  having  a  side 

—1  —2  _8 

of  length  2,2  ,*..,2  respectively, 

—8 

SMALLER  SIZES:  gives  the  number  of  elements  smaller  than  2  ,  and 

MAX. AD J. RATIO:  gives  (in  power  of  2)  the  maximum  of  ratios  of  sizes 

of  two  adjacent  elements  in  the  2-D  domain. 

The  next  line  gives  the  sum,  maximum,  mean,  and  threshold  value  for  the 
error  indicators  of  the  elements  in  the  2-D  domain. 

The  following  three  lines  give  statistical  information  on  the  error 
indicators.  The  value  in  parenthesis  gives  the  total,  and  the  following  8 
numbers  give  the  number  of  elements  in  8  standard  deviations  from  the  mean 
or  maximum  predicted  value. 

The  line  which  gives  the  TIME  breakdown  is  discussed  in  Section  4.4. 
After  the  storage  requirements  are  printed,  a  list  of  the  nodal  points 
of  the  mesh  in  the  2-D  domain  are  given.  Included  are  the  index  of  the  point 
as  determined  from  the  numbering  order  described  earlier,  as  well  as  the 
point's  local  coordinates,  global  coordinates,  computed  solution  values, 
and  the  error  of  the  computed  solution  (if  the  exact  solution  is  known  and 
defined  in  ZMPTRU) .  Only  regular  points  are  listed. 

After  the  nodal  points  are  listed,  the  elements  are  listed  in  order 
of  decreasing  error  indicators.  The  indices  of  the  elements  are  determined 
through  the  same  numbering  system  as  the  nodal  points  and  this  numbering 
order  was  described  earlier.  Also  included  are 
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— H+l 

H  -  indicating  the  size  of  the  element  (2  )  » 

R  -  indicates  the  regularity  of  the  cornerpoints  of  the  element,  e.g. 

R-3  corresponds  to  all  regular  cornerpoints  (see  Section  4.2). 
LOCAL  COORD  -  the  local  coordinates  of  the  center  of  the  element, 
ERROR  IND.  -  the  error  indicator  for  the  element,  and 
PREV.ERR.IND.  -  the  previous  error  indicator.  See  Appendix  A.3  for 
further  details  on  the  previous  error  indicator. 

4-4  The  TIME:  breakdown 

FEARS  uses  substructured  solving,  obtaining  solutions  first  for  those 
nodes  on  the  1-D  and  0-D  domains  and  then  backsolving  to  obtain  solutions 
for  the  nodes  in  the  interior  of  each  2-D  domain.  In  order  to  better 
understand  the  TIME  breakdown,  write  the  assembled  global  stiffness  matrix 
in  the  form: 


where  X.  correspond  to  the  unknowns  in  the  interior  of  2-D  domain  i,  and  X. 

i  BD 

the  unknowns  for  the  0-D  and  1-D  domains. 

With  the  2-D  domain  print  command,  the  line  starting  with  TIME  lists 
ASM  &  DEC:  the  time  for  the  assembly  and  LU  decomposition  of  A^ 

BCK:  the  time  for  the  backsubstitution  for  the  unknowns  in  2-D  domain  i. 

THR:  the  time  for  computing  the  error  indicators,  and 

TOTAL:  the  sum  of  the  above  three  times. 
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The  REPORT  print  lists 
EXECUTION  TIME: 

SUBDIVISION:  the  time  of  subdivision  for  the  full  domain, 

2-D  MATRIX  SOLUTION:  the  sum  of  the  ASM  &  DEC  times  described  above  for 
each  2-D  domain, 

BDRY  MATRIX  SOLUTION:  the  time  for  solving  C  •  X^  -  YBD  where  C  ,  YflD 

are  C  ,  modified  by  the  2-D  partial  decompositions. 

2-D  ERROR  CALCULATION:  The  sum  of  the  THR  times. 

At  the  end  of  an  AUTO  command  the  ***AUTO  TIME  gives  the  sum  of  the  4 
execution  times  plus  overhead. 


4.5  The  REPORT  printout 

The  REPORT  has  the  format 
*****  FULL  DOMAIN  ***** 

NUMBER  OF  POINTS:  NP  NUMBER  OF  ELEMENTS:  NE 

ENERGY  NORM:  NM  ENERGY:  ENG 

ERROR  ESTIMATOR:  EST  RELATIVE  ERROR:  REL 

MAX. ERROR  INDICATOR:  MAX  BY  2-D  INDEX _  THRESHOLD:  TH 

APPROXIMATE  NUMBER  OF  ELEMENTS  TO  BE  SUBDIVIDED: 

2-D  NO.  OF  ELEMENTS 


TOTAL  T 

STORAGE  SIZES;  MAX, CORE- _ ,  TOTAL- _ ,  NO  RECORDS-_ 

BDRY  MATRIX  -  ___ 

EXECUTION  TIME:  ,  SUBDIVISION: _ 2-D  MATRIX  SOLUTION: 


BRDY  MATRIX  SOLUTION: _ ,  2-D  ERROR  CALCULATION: _ 

Here 

NUMBER  OF  POINTS:  gives  the  number  of  nodes  in  the  0,  1  and  2-D  domains 

have  at  least  one  degree  of  freedom, 

NUMBER  OF  ELEMENTS :  gives  the  total  number  of  elements  in  the  full  domain, 
ENERGY  NORM:  NM  -  VENG  where 

ENERGY:  ENG 

ERROR  ESTIMATOR:  gives  the  error  estimate  for  the  full  domain, 

RELATIVE  ERROR:  REL  -  NM/EST. 

The  next  line  gives  the  maximum  error  indicator,  the  2-D  domain  it  is  in, 
and  the  threshold  value  for  automatic  refinement. 

Next  is  a  list  of  the  number  of  elements  in  each  2-D  domain  which  will 
be  subdivided  if  the  threshold  value  is  used  for  subdivision. 

Storage  size  information  is  then  given  and  finally  a  breakdown  of  the 
times  as  discussed  in  the  previous  section. 

4.6  The  OUTPUT  printout 

In  this  section  we  describe  the  data  either  printed  or  file  written, 
with  either  the  OUTPUT  command  or  by  previously  setting  IPR(8)  in  the  DEBUG 
command  or  at  the  program  start.  The  data  printed  will  be  a  list  with  the 
headings. 

2-D  ELT  H  LOCAL  COORD  GLOBAL  COORD  ERR. IND  0UTP(1)  .  .  .  0UTP(5) 

The  column  headed  with 

2-D  -  gives  the  index  of  the  2-D  domain  that  the  element  is  in, 

ELT  -  gives  the  index  of  the  element, 

—H+l 

H  -  indicates  the  element  size  (-2  ), 

LOCAL  COORD  GLOBAL  COORD  -  gives  the  local  and  global  coordinates  at  the 
center  of  the  element. 
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ERRIND  -  gives  the  error  indicator  of  the  element, 
0UTP(1) 


0UTP(2) 

0UTP(3) 

0UTP(4) 

0UTP(5) 


^  will  give  the  stresses  o^,  a^t  o^,  and  solutions  u^,  U2  at 
center  of  the  element,  if  S  was  given  for  the  elasticity  problem. 


These  last  5  values  are  determined  by  the  output  matrix  S  given  in  the  bilinear 
matrix  input.  The  5  values  printed  are  determined  by 


ru 

1 

§ 

3 

h* 

32 

... 

LuJ 

_0UTP(5)J 

6x1 

• 

5x1 

The  appropriate  s  for  elasticity  problems  is  given  in  the  Appendix. 


4.7  DUMP-File  output 

A  sequential  binary  DUMP-file  is  generated  by  the  DUMP  command.  This 
file  can  be  used  as  input  for  FEARS  (see  Sections  3.13  and  3.20),  or  as  input 
for  various  postprocessors.  The  order  of  the  binary  records  on  the  file  are 
as  follows: 

1.  Summary  record  of  the  problem  (16  words) 

2.  Last  long  path  history  record  (40  words) 

3.  Function  parameter  record  (N+l  words  where  N  is  the  number  of 
parameters) 

4.  0-D  summary/records  (8  words /records,  Nq  records) 

5.  1-D  summary  records  (18  words/records,  records) 

6.  2-D  Summary  records  (25  words/records,  ^  records) 

7.  For  each  1-D  domain  with  non-fixed  boundary  condition: 

7,1  l-D  data  record  (8k  words,  where  k  is  the  number  of  points 
on  the  1-D) 


L 
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7.2  If  Y»e  matrices  were  defined  for  the  1-D  then  data 


record  (13  words) 

8.  For  each  2-D  domain 

8.1  2-D  data  record  (4k  words,  where  k  is  the  number  of  points 
and  elements  in  the  2-D) 

8.2  Bilinear  form  data  record  (125  words) 

Details  on  the  fields  of  records  can  be  found  in  the  Program  Documentation, 


CHAPTER  V.  STARTING  THE  FEARS  PROGRAM 


This  chapter  describes  the  UNIVAC  1100  EXEC  control  cards  necessary 
to  run  the  FEARS  program.  There  are  two  elements  and  one  file  which  are 
of  concern: 

Absolute  element:  CKM*FA. A 

Relocatable  file:  CKM*RE. 

Map  element:  CKM*SE.MAP 

The  names  of  the  above  elements  (file)  may  change  by  the  actual  installation. 
The  absolute  element  contains  the  absolute  program  with  dummy,  zero  valued, 
function  routine  for  E(X,Y) ,D(S,Y)  which  can  be  run  if  that  routine  is 
satisfactory  by 

@XQT,F  CKM*FA.A 

Otherwise,  section  5.2  describes  how  to  set  up  an  absolute  program  with  the 
help  of  the  Relocatable  file  and  Map  element  to  incorporate  the  user's 
supplied  function  routine. 

Section  5.1  explains  the  necessary  file  assignments  to  be  included  prior 
to  the  run  of  the  program.  Section  5.3  describes  the  necessary  initial 
inputs  for  FEARS  before  the  geometry,  bilinear  matrices  and  command  inputs 
are  given  as  described  in  chapters  1,  2  and  3,  respectively. 

5.1  File  Assignment 

FEARS  uses  six  temporary  files  (11  to  16)  with  all  runs  which  should  be 
assigned  by 

@ASG,T  11. 

@ASG,T  12. 

@ASG,T  13. 

@ASG,T  14. 

@ASG,T  15..///512 

@ASG,T  16. 

Four  other  files  may  be  used  depending  on  the  actual  run  and  commands  given 
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for  the  run.  These  files  should  be  catalogued  files  with  user's  given  names. 

These  names  should  be  associated  with  the  fortran  unit  number  by  @USE: 

@ASG,  A  file  name. 

@USE  n. ,  file  name. 

where  n  is  the  FORTRAN  unit  number.  The  four  files  are  as  follows: 

(i)  Output  file  -  If  the  OUTPUT  command  is  used,  see  3.11,  or  IPR(8)  was 

given  as  2  or  3  in  either  the  DEBUG  command  (see  3.10) 
or  initially  (see  5.3),  the  file  should  be  assigned  with 
n*17  .  At  the  termination  of  the  run,  the  file  will 
contain  the  stresses  as  described  4.6. 

(ii)  Mesh  file  -  This  file  is  generated  as  output  file  if  IPR(l)  ■  -n  as 

given  for  the  DEBUG  command  (3.10),  it  is  used  as  input 
file  for  MESH  command  (3.19;  It  is  recommended  to  use 
n  ■  18. 

(iii)  Dump  file  -  This  file  is  generated  by  the  DUMP  command  (3.12).  Since 

present  and  future  postprocessors  use  this  file  as  input 
file,  most  of  the  use  of  FEARS  is  anticipated  to  use  this 
file.  Recommended  n  ■  20. 

(iv)  Reset  file  -  This  input  file  is  assumed  to  be  generated  by  a  DUMP 

command  in  previous  run  of  FEARS.  The  file  is  used  either 
initially  (5.3)  or  by  the  RESET  command  (3.13).  Recommended 
n  ■  21. 


r 


L1 


l: 


5.2  Preparation  for  Execution 

When  the  user  has  to  incorporate  his/her  function  routine  (see  2.7),  then 
a  new  absolute  program  of  FEARS  must  be  generated.  The  recommended  steps  to 
be  taken  are  as  follows: 

(i)  Write  the  function  routine  according  to  the  spef ideation  given  in  2.7. 
Compile  the  symbolics  using  the  ASCII  FTN  compiler  to  generate  the 
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relocatable  element.  For  further  on,  it  is  assumed  that  F.FXY  is 
element  name  of  this  relocatable  element. 

(ii)  Using  the  editor, 

@ED  CKM*SE.MAP,M 
change  the  Map  element  into  M 
user  *N 

computer  @MAP,IN  ,CKM*FA.A 

user  *C  /CKM*FA/F/ 

computer  @MAP,IN  ,F.A 

user  *n  2 

computer  IN  CKM*RE .  FUNC 

user  *C  /CKM*RE . FUNC/F . FXY/ 

computer  IN  F.FXY 

user  *@ 

computer  editor  signs  off 

user  @ADD  M 

computer  MAP  ...  . 

END  MAP.  ERROR  -  0  ... 

The  above  sequence  places  the  new  absolute  element  F.A  in  the  users 
file. 

(iii)  After  the  appropriate  assign  statements,  see  5.1,  the  run  can  be 
initiated  by 


The  geometry  and  matrix  Inputs  can  be  prepared  separately  in  file 

elements,  e.g.  in  INP.G1  and  INP.MX1  ,  respectively.  In  this  case, 

the  user  simply  applies  the  UNIVAC  commands 

@ADD  INP.G1 
@ADD  INP.MX1 

in  the  above  input  stream  as  input  for  the  geometry  and  bilinear,  error 
and  output  matrices.  This  technique  is  especially  useful  when  FEARS  is 
used  interactively. 


5 . 3  Execution  of  FEARS 


After  the  program  execution  has  been  initiated  by  @XQT,F,  the  computer 

will  respond 

*****  FEARS  ***** 

2-D  FINITE  ELEMENT  PROGRAM 
UNIVERSITY  OF  MARYLAND.  1981. 

MAXIMUM  ALLOWED  SIZES 
NUMBER  OF  0-D,  1-D,  2-D:  34  49  16 
NUMBER  OF  POINTS  ON  A  1-D:  31 
NUMBER  OF  POINTS  AND  ELEMENTS  IN  ONE  2-D:  482 
MATRIX  STACK  SIZE  FOR  ONE  2-D:  8800 

BOUNDARY  MATRIX  SIZE:  14000 

SHORT  PATH  ERROR-FACTOR  -  .55000 

PRINT  CONTROL  INTEGERS  (8): 

IPR(l)  -  PRINTS  DURING  SHORT  PATH,  SUBDIV. 

NEG.  K  -  RECORD  SUBD.  ELEMENTS  ON  FILE  K 
IPR(2)  -  PRINTS  DURING  ERROR  CALCULATION 

IPR(3)  ~  PRINTS  DURING  ASSEMBLY 

IPR(4)  -  PRINTS  DURING  DECOMPOSITION 

IPR(5)  -  PRINTS  DURING  MATRIX  SOLUTION 

IPR(6)  -  PRINTS  DURING  BACKSUBSTITUTION 

IPR(7)  -  REPRINTS  INPUT 

IPR(8)  -  AUTOMATIC  ELEMENTAL  OUTPUT  CONTROL 
IPR(J) ,J«1, . . . ,8: 

The  user  now  should  input  the  8  integer  values  IPR(l)  -  IPR(8)  as  discussed 
in  the  DEBUG  command,  e.g. 

0,0, 0,0, 0,0, 1,0 

Next  the  computer  will  respond 

PROBLEM  INPUT  DATA  IS  ON 
FILE  NUMBER  : 

The  user  should  input  zero: 

0 

if  this  is  a  new  problem,  or  the  Reset  file  number,  e.g. 

21 

on  which  all  the  data  at  some  stage  of  a  problem  was  saved  by  the  DUMP  command. 
In  this  case,  the  computer  will  respond 
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****  COMMAND 


£ 


6 


ff 


» 


H 


and  the  prograa  will  be  in  the  command  mode. 

Otherwise  (in  case  of  a  new  problem) ,  the  computer  will  respond 
FUNCTION  PARAMETERS:  N,P1, . . . ,PN(N>0) : 
and  the  user  should  input  an  integer  N>0  ,  and  N  parameters  for  use  in  the 
user  defined  functions  routine  (see  section  2.7).  If  no  parameters  are 
required,  the  dummy  values  1,0,  may  be  input. 

The  computer  will  then  respond 

ID  -  NUMBER  OF  THE  PROBLEM: 

and  any  integer  response  from  1  to  999999  will  suffice. 

The  next  computer  response  is 

P-NORM  FOR  THE  FULL  DOMAIN: 

The  P  requested  here  is  the  same  one  as  in  (2.1).  If  an  L_  norm  is  of 

Zp 

interest  in  measuring  the  errors  then  the  value  P>.5  should  be  input.  An 
Lm  norm  will  be  used  if  0.  is  input.  Usually  P*1  which  corresponds  to 
the  standard  L2  energy  norm. 

The  computer  will  then  respond 

PROBLEM  ID:  _ 

DATE:  _ 

P-NORM:  _ 

GEOMETRY :  _ 

At  this  point  the  geometry  should  be  input  in  the  format  described  in  chapter 
1,  i.e.,  starting  with  number  of  0-D  domains  and  ending  with  the  last  2-D 
domain  description  line. 

If  no  errors  are  detected  the  computer  will  respond 
GEOMETRY  ACCEPTED 

BILINEAR,  ERROR,  AND  OUTPUT  MATRICES: 
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At  this  point,  these  matrices  should  be  input  in  the  format  described  in 
Chapter  2. 

The  computer  should  then  respond 
MATRICES  ACCEPTED 
INITIAL  SUBDIVISION,  1  OR  2: 

The  value  1  will  cause  each  2-D  domain  to  have  4  elements  and  the  value 
2  will  cause  each  2-D  domain  to  have  16  elements  for  the  initial  mesh. 
Finally  the  computer  will  respond 

INITIAL  SUBDIVISION  (  )  FOR  1-D  PERFORMED 
INITIAL  SUBDIVISION  (  )  FOR  2-D  PERFORMED. 

The  program  will  then  obtain  an  initial  full  solution  path  on  the  entire 
domain,  print  a  REPORT  and  respond 
****  COMMAND 

signifying  that  it  is  ready  for  user  commands. 

Once  in  conmand  mode,  further  actions  are  governed  by  the  user's 
commands  as  described  in  Chapter  3. 
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APPENDIX  A 


A.l  Sample  Geometry 

In  this  appendix  sample  domains  and  their  FEARS  geometry  input  are 
presented. 

Example  1:  Quarter  Ring  -  1  -  subdomain 


Figure  A.l.  Quarter  Ring 
Trie  geometry  input  for  Figure  A.l  is: 

4 

1,0. , .5,1,0. ,0. 

2,0. ,1. ,1,0. ,0. 

3. . 5.0. .2.0. .0. 

4.1. . .0. .2.0. .0. 

4 

1,1, 3, 0,2. 

2, 2, 4, 0,1. 

3, 1,2, 1,0. 

4, 3, 4, 2,0. 

1 

1,1, 2. 3,4 

Note  that  the  boundary  condition  on  line  index  _1  is  specified  as  free, 
force  present  (Indicated  by  arrows)  is  specified  in  the  bilinear  matrix 
input. 


The 
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Example  2:  Unit  Square  -  4  -  sub  domains 


The  geometry  input  for  Figure  A. 2  is: 


9 


1,0. ,0. ,3,0. ,0. 

4,5, 8, 0,0 

2,0*  y .5,3,0. ,0. 

5, 1,2, 3,0 

3,0. ,1* ,3,0. ,0. 

6, 4, 5, 0,0, 

A, .5,0. ,0^0. ,0. 

7, 7, 8,3,0 

5, *5, .5,0,0. ,0. 

8, 2, 3, 3,0 

6, .5,1* ,0,0. ,0. 

9, 5, 6, 0,0 

7,1. ,0# ,3,1. ,0. 

10,8,9,3,0 

8,1. ,.5,3,1. ,0. 

11,3,6,0,0 

9,1. ,1* ,3,1. ,0. 

12,6,9,0,0 

12 

4 

1,1, 4, 0,0. 

1,1, 2,4*5 

2, 4, 7,0,0. 

2, 2, 3, 5, 6 

3, 2, 5, 0,0. 

3. 4. 5. 7. 8 

4. 5. 6. 8. 9 

Example  3 


Figure  A. 3. 

Disk  with  hydrostatic  force 

5, 5, 1,0,0. 

,0. 

6, 6, 2, 0,0. 

). 

7,3, 7,0,0. 

5. 

8, A, 8, 0,0. 

,0. 

9, 8, 5, 0,1. 

,0. 

10,5,6,0,1. 

). 

11,6,7,0,1. 

). 

12,7,8,0,1. 

►0. 

5 

1.1. 2.3. 4 

2. 5. 1.8. 4 

3. 5. 6. 1.2 

4, 2, 6, 3, 7 

5. 8.4. 7. 3 

Notice  the  boundary  conditions  at  points  5  and  7.  For  this  problem, 
additional  boundary  conditions  are  necessary  to  ensure  uniqueness  of  a 
solution  by  eliminating  rotations  and  translations  of  the  solution.  Poini 
7  was  picked  arbitrarily  as  the  stationary  point  (u^=U2=0) .  To  prevent 
rotations  the  vertical  displacement,  ,  was  set  to  zero  at  point  5.  It 
would  be  equivalent  to  fixing  the  vertical  displacement  at  either  point  1 
or  point  3  as  well. 


Notice  that  points  2  and  3  have  the  same  coordinates  and  thus  lines 


2  and  3  lie  on  top  of  each  other.  This  represents  the  two  edges  of  a 
crack. 


A.  2  Bilinear  and  Error  Matrices 


A. 2.1.  The  Matrix  A  Arising  From  Elasticity 

The  principle  of  virtual  work  for  the  elasticity  equations  yields 


an  integral  of  the  form  J  (6e)AC£  ,  where 

D 


■m 

3ui 

XX 

3x 

m 

3u2 

yy 

3y 

3u. 

xy 

1 

ay 

and  C 


is  the  3x3  stress-strain  matirx  such  that 
"o 

“I 

e  .  For  the  plane  strain  assumption  the 


yy 


xy 


matrix  C  is  given  by 


E 


(l+v)(l-v) 


1-v  v  0 

v 

0  0 


1-v  0 

l-2v 


and  for  the  plane  stress  assumption 

„  _  E 


1-v 


1 

v 

0 


v  0 


1  0 
«  1-v 


where  E  is  Young's  raolulus  and  v  is  the  Poisson  ratio. 
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The  corresponding  integrand  in  the  FEARS  formulation  is 


*3Z 


Since 


au 

az 


au. 


ax 

3U„ 


ax 

3D, 


av 

au. 


ay 


write 


3U 


az 


,  where 


.av 


1 

o 

0 


0 

0 

1 


0  0 
0  1 
1  0 


Also,  Sc  =  Djj^  •  Therefore, 


T 

(S_e)TC  e  *[H]  °T  C  D[|z]  ,  from  which  we  have 


A  -  D  C  D  . 

For  the  plane  strain  assumption 

1-v  0  0  v 


(1+v)  (1-  2v) 


0 


lz2v  1^2v 
0  2  2  u 


-v 


and  for  the  plane  stress  assumption 


1  0  0  v 

«  l-2v  l-2v  A 


we  may 


the  relationship 


A. 2.2  The  Error  Matrix 


For  elasticity  problems  the  matrices  (C)j  described  in  Section  2.3 
are  zero,  and  so,  from  (2.1)  the  error  is  approximated  in  the  norm 


lllulil 


2p 


A  UBKffl} 


j 


-hr* 


To  obtain  error  estimates  for  the  energy  norm,  simply 


take  (Ag)  *  (A)^  . 
Suppose,  however,  that  our  interest  is  concentrated  on  the  error  in 


one  of  the  stresses,  say  o 


xx 


Since  a 


2  T 
a  **  a 
xx 


[l  0  o]  a 


XX 

_t  r 


[i  0  o]2  , 


10  0 
0  0  0 
0  0  0 


T_T 

£  C 


1 

0 

0 


0 

0 

0 


0 

0 

0 


C  e 


m 


T  T 
D  C 


1 

0 

0 


0 

0 

0 


0 

0 

0 


C  D 


[f]  . 


where  C  is  the  3x3  stress  strain  matrix,  and  D  is  defined  by  (A.l). 
Let  Q  ■ 


1 
0 
0 

we  should  take 


0  0 
0  0 
0  0 


Then,  if  the  error  a  is  our  main  concern, 

xx 


T  T 

Ag  -  D1  CA  Q  C  D  . 

For  example,  using  the  plane  strain  assumption  with  an  error  in  , 


*E  “  [(1+v)  (1-v)] 


(1  -vr 
0 
0 

v(l-v) 


0 

0 

0 

0 


0 

0 

0 

0 


v(l-  v) 
0 
0 
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For  the  error  in  a ^  take  Q  in  (A.  2)  to  be 


for  o  ,  take  Q  to  be 
xy 


Q  -  0  0  0 

0  10 


0  0  0  I ,  and 


0  0  0 

0  0  0 

0  0  1. 


A. 2. 3  The  Output  Matrix  S 

For  the  OUTPUT  printout  (see  Sections  3.11  and  4.6)  the  vector  printed 

is  obtained  from  the  multiplication 

S  3U  where  3U  and  U  are 
I  3Z  I  3Z 


evaluated  at  the  center  of  each  element.  In  order  for  a  printout  of 
°xx*  cyy’  axy*  ul’  u2  the  matrix  s  should  be 


1-v  1— v 

2  2 

0  0 


under  the  plane  stress  assumption,  and 


(1+v) (l-2v) 


0  1-v  0 


l-2v  l-2v 


under  the  plane  strain  assumption 


Example  1 


A-Matrix:  Laplacian 
Error  norm:  H^-  seminorm. 


Packages:  One  package  for 
3u.  3u2 

Output  Matrix:  -gj-,  -gj-. 


Line  Integrations:  None 


V 


Example  2:  A  Matrix:  Plane  Strain  Elasticity 
Error  Norm:  Energy 

Packages:  Two  packages  for  4  2-D  domains. 

Package  1:  2-D  domains  1  and  2,  v  *  0.0,  E  *  1.0 
Package  2:  2-D  Domains  3  and  4,  v  *  0.3,  E  ■  1.0. 
Output  Matrix:  o^,  a  ,  o^,  u2,  u2* 

Line  Integrations:  Local  normal  force  on  line  4. 


2 

1,  2,  1,  2,  0.,  1,  0,  0,  0,  0 

1.,  0. ,  0.,  *3,  0.,  0.,  .3,  0.,  0.,  .3,  0.,  0.,  .3,  0.,  0.,  1. 
I* »  0* »  0* »  0.,  0.,  .5,  0.,  0.,  .5,  0.,  0.,  .5,  0.,  0.,  1. 

1. ,  1. ,  1. ,  1. 

1.,  0.,  0.,  0.,  0 « ,  0. 

0. ,  0.,  0.,  1.,  0.,  0. 

0.,  .3,  .3,  0.,  0.,  0. 

0.,  0.,  0.,  0.,  1.,  0. 

0. ,  0.,  0.,  0.,  0.,  1. 

2,  2,  3,  4 
1,  0,  0,  0,  0 

1.3461538,  0.,  0.,  .38461538,0. ,  .57692308,  .38461538,  0. 

0.,  .38461538,  .57692308,  0.,  .38461538,  0.,  0.,  1.3461538 

1.3461538,  0.,  0.,  .38461538,  0.,  .57692308,  .38461538.0. 

0.,  .38461538,  .57692308,  0.,  .38461538,  0.,  0.,  1.3461538 

1. ,  1. ,  1. ,  1. 

1.3461538,  0.,  0.,  .57692308,  0.,  0. 

.57692308,  0.,  0.,  1.3461538,  0.,  0. 

0.,  .38461538,  .38461538,  0.,  0.,  0. 

0.,  0.,  0.,  0.,  1*,  o. 
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A. 3  Computation  of  the  Threshold  Value 


When  using  the  AUTO  command  or  the  ITERATE  command  (the  ITERATE  command 
causes  a  sequence  of  AUTO  commands  to  be  performed) ,  it  is  Important  to 
understand  how  the  threshold  value  is  computed.  Recall  that  with  the  AUTO 
command,  all  elements  having  error  indicators  above  a  certain  threshold 
are  subdivided.  Because  of  the  remarks  made  in  Section  3.1,  we  seek  a  mesh 
in  which  all  error  indicators  are  nearly  equal. 

The  most  naive  way  to  equalize  the  Indicators  would  be  to  simply  sub¬ 
divide  the  element  with  the  largest  error  indicator  before  each  solution 
path.  Suppose  that  the  elements  in  the  initial  mesh  are  labeled 

1,  2, . . .  ,N 

in  order  of  decreasing  error  Indicators.  Following  our  "naive"  strategy, 
we  sould  subdivide  element  1  into  four  subelements — 1.1,  1.2,  1.3,  and 
1.4.  After  resolving  and  recomputing  the  error  indicator  our  new  list  might 
be 

2,  3, . . . J-,  ,  1*1,  *1  ^ 1 1, . . . , Jj*  1.2,  J 1, . . .  ,N  . 

Next,  element  2  gets  subdivided  into  2.1,  2.2,  2.3,  2.4,  etc.  Suppose  element 

1.1  has  the  largest  indicator  of  the  first  4J^  subdivided  elements.  That  is, 

♦ 

E(l.  1)  E(i.  j)  for  1  ■  2, . . .  ,«J^ *  and  j  ■  1 , . . . , 4  . 

Then  our  strategy  would  have  subdivided  the  first  elements  one  at  a  time, 

obtaining  a  new  solution  between  each  subdivision.  Clearly,  it  would  be  more 
efficient  to  subdivide  all  elements  at  once  before  obtaining  a  new 

solution.  The  cufoff  at  could  easily  be  determined  if  we  know  what  the 

maximum  error  indicator  would  be  for  the  next  level  of  subdivided  elements 
(in  this  case  E(l.l)).  The  following  procedure  is  used  to  approximate  this 
cutoff  value. 
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Suppose  that  a  father  element  F  is  to  be  subdivided  into  its  four 
sons — F.l,  F.2,  ^.3,  and  F.4.  (see  Figure  A. 5) 


F.2 

F.4 

F.l 

F.3 

Figure  A. 5. 

The  error  indicator  for  the  father  element,  E(F)  is  then  saved  as  the 
previous  error  indicators  for  the  four  sons.  It  is  then  assumed  that  the 
errors  arising  from  subdividing  the  sons  F.l,  F.2,  F.3,  and  F.4  will 
decrease  by  the  same  ratio  as  from  the  previous  subdivision.  That  is,  we 
use  the  recipe:  PREDICTION  -  x  PRESENT.  Thus,  the  predicted  error 

indicator  for  each  of  the  four  sons  is  computed  by  the  formula 

P(F.i)  -  E(F. i)2/E(F)  ,  for  i  -  1,  2,  3,  4  . 

In  order  to  ensure  that  the  predicted  error  is  smaller  then  the  present 
error  we  add  the  condition— 

P(F.i)  -  min/P(F.i),f  (.9)2pE(F.i)  if  p>  .5^ 

\  \(.9)E(F.i)  if  p  -  0  J)  , 

where  p  is  the  p-norrn  for  the  2-D  domain  as  input  in  the  (W  )  line  of 

c  i 

bilinear  and  error  matrix  input. 

If  an  element  is  not  a  son  of  some  father  element,  that  is,  the  element 
also  belongs  to  the  initial  mesh,  then  no  previous  error  indicator  is 
available.  In  this  case,  we  make  a  prediction  for  the  error  indicators 
upon  subdivision  of  this  element  F  ,  by  the  formula 

P(F)  -  r (/2?2)2pE(F)  if  p>  .5 

\</272)E(F)  if  p  -  0  . 

The  threshold  T  is  then  calculated  by  the  formula 
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T  -  max  P(t)  , 

T 

as  T  ranges  over  all  elements  In  the  mesh 


i 


i  p 

m 
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APPENDIX  B 


FEARS  Element  Postprocessor 

The  purpose  of  the  Postprocessor  is  to  allow  flexible  computations  of 
various  functionals  on  the  approximate  solution.  The  functionals  can  have 
various  structures  and  can  be  used  for  the  effective  computation  of  the  stress 
intensity  factors  in  fracture  mechanics  etc.  At  the  present  the  functionals 
are  on  the  level  of  2-D  domains,  but  will  be  extended  in  the  future  to  deal 
with  functionals  on  the  1-D  domains  also. 

A  FEARS  Postprocessor  frame  has  been  set  up  in  the  file  CKM*ELTOUT. 

Both  symbolic  and  relocatable  elements  for  an  Element  Postprocessor  has 
been  established  in  the  file.  To  produce  an  executable  absolute  element, 
the  user  must  write  some  subroutines  and  combine  it  with  the  above  file. 

The  FEARS  Element  Postprocessor  reads  the  RESET  file,  generated 
by  the  FEARS  program  using  the  DUMP  command.  The  Postprocessor  requires  a 
user's  written  subroutine  (with  5  entries).  The  purpose  of  the  Postprocessor 
frame  is  to  relieve  the  user  of  the  cumbersome  task  of  setting  up  the  data 
structure  from  the  RESET  file  and  setting  the  individual  element  informations 
needed  for  his/her  calculations.  The  frame  allows  it  to  process  all  elements, 
or  only  those  which  are  in  certain  2-D  domains  designated  as  active  2-D 
domains  by  the  user. 

Frame  algorithm: 

The  main  program  of  the  Postprocessor  frame  has  the  following  algorithm: 

1.  Read  the  RESET  file  and  set  up  Internal  data  structure. 

2.  Clear  activity  tags  MDlTMP(j)  and  MD2TMP(k)  of  all  1-D  and 
2-D  domains,  repectively. 

3.  Call  ELTO(MSARY)  calls  user's  initialization  routine. 

4.  DO  1*1  to  MSP2  (number  of  2-D). 

5.  IF  MD2TMP(I)-0  THEN  GOTO  Next  I 


6.  Set  up  core  storage  for  2-D  index  I. 

7.  CALL  ELT1(MD2ARY) . 

call  user's  2-D  initialization  routine 

8.  DO  J*1,M2SXE  (number  of  elements  in  the  2-D) 

9.  CALL  ELTX(K,N,X,C) 

call  user's  processing  routine 

10.  Next  J 

11.  CALL  ELT1E (MD2ARY) 

call  user's  2-D  summary  routine 

12.  Next  I 

13.  CALL  ELTO(HSARY) 

call  user's  final  summary  routine 

14.  STOP 


User’s  supplied  subroutine (s) : 


The  user  must  supply  the  following  five  subroutines,  or  subroutine(s) 
with  the  appropriate  entries.  Note  that  the  procedures,  ZMS,  ZMDO,  ZMD1 
and  ZMD2  may  be  used  in  these  routines.  These  procedures  define  data  areas, 
summary  records  of  the  problem.  The  user  may  also  use  procedures  ZP2T, 
ZP20S,  ZP21S  and  ZP2B,  which  give  the  data  area  (restored)  for  a  given  2-D 
domain.  For  proper  use  of  these  data,  the  user  should  consult  FEARS  program 
documentation . 

1.  SUBROUTINE  ELTO (MS ARY) 


INCLUDE  ZMS 


INCLUDE  ZND1 
INCLUDE  ZMD2 


This  routine  is  called  as  an  initialization.  The  user  must  set  the 


active  2-D  domains  by  setting  MD2TMP(j)  to  non-zero  with  j  being  the 
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active  2-D  domain  indices.  These  indices  may  be  obtained  by  user's 
defined  input,  or  in  case  of  all  2-D's,  through  a  simple  loop  on  j=l,MSP2, 
where  MSP2  is  defined  in  procedure  ZMS  (summary  array  of  the  full  domain) 
and  MD2ARY  is  defined  in  procedure  ZMD2  (summary  records  of  2-D  domains) . 
user  may  also  mark  the  1-D  domains,  if  needed  in  later  calculations.  This 
can  be  done  by  storing  Integer  values  in  MDlTMP(k)  where  k  is  the  index 
of  the  1-D,  and  MD1ARY  1-D  summary  records  are  defined  in  procedure  ZMD1.  On 
subsequent  processing,  these  integer  values  may  be  tested. 

The  user  may  define  further  storage  areas  needed  for  subsequent 
calculations,  also  perform  any  input  needed,  outputs  such  as  headings,  etc. 

2.  SUBROUTINE  ELT1(MA2) 

DIMENSION  MA2(26) 


This  routine  is  called  in  the  beginning  of  all  active  2-D  domains  where 
the  array  MA2  is  the  summary  array  of  the  2-D  domain  with  MA2(1)  containing 
the  index  of  the  2-D,  MA2(26)  the  activity  non-zero  tag  set  by  the  user. 

3.  SUBROUTINE  ELTX(K,N,X,C) 

DIMENSION  X(2) ,C(2,4) 


This  routine  is  called  for  each  element  in  the  active  2-D  domains  where 
the  user  must  perform  his  own  calculations.  The  arguments  are  input  arguments 
as  follows: 

K  ■  index  of  the  element 

N  -  level  number  of  the  element  in  the  transfered  unit  square, 
side  length  h«2**(-N-l) 

X(l) ,X(2)  ■  local  (unit-square)  coordinate  values  of  the  middle 
point  of  the  element 


C(l,j),C(2,j)  ■  solution  values  at  the  j'th  cornerpoint  of  the 


element  where  j  is  defined  in  order: 

2 - 4 

•  • 

•  • 

1 - 3 

Furthermore,  the  following  function/subroutines  are  also  available  for  use: 
H  -  FUTH(N) 

side  length  of  the  element 
CALL  ZP2XY (X,Z) 

gives  the  global  coordinates  Z(1),Z(2) 
corresponding  to  the  local  coordinates  X(1),X(2) 

CALL  ZP2TRX(X,T,D) 

gives  transformation  matrix  T  and  its  determinant  D  at  the 
local  coordinates  X(l) ,X(2) : 

T(l)  -  d(Z(l))/d(X(l)) 

T(2)  -  d(Z(l))/d(X(2)) 

T(3)  -  d(Z(2))/d(X(l)) 

T(4)  -  d(z(2))/d(X(2)) 

4.  SUBROUTINE  ELT1E (MA2) 

DIMENSION  MA2(26) 


This  routine  is  called  after  all  elements  of  an  active  2-D  domains 
have  been  processed  by  the  ELTX  routine.  This  entry  allows  the  user  to 
print  any  summary  for  the  2-D  domain,  MA2  array  is  the  same  as  for  the 
subroutine  ELT1. 

5.  SUBROUTINE  ELTO(MSARY) 

DIMENSION  MS ARY (16) 
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0 


0 
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END 
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